#include "BVP.h"
#include <filesystem>
#include <fstream>
#include <jsoncpp/json/json.h>

#define SECTOR "-----------------------------------"

using namespace std;
namespace fs = std::filesystem;

double f(const Point2d &p)
{
    double x = p[0], y = p[1];
    return (sin(x) - cos(x) * cos(x) - 1) * exp(y + sin(x));
}

double u(const Point2d &p)
{
    double x = p[0], y = p[1];
    return exp(y + sin(x));
}

double u_x(const Point2d &p)
{
    double x = p[0], y = p[1];
    return cos(x) * u(p);
}

double u_y(const Point2d &p)
{
    return u(p);
}

double g_partial(const Point2d &p, const BVP2d::Disk &D)
{
    Point2d c;
    double r;
    if (D.radius > 0)
    {
        c = D.center;
        r = D.radius;
    }
    Point2d partial{u_x(p), u_y(p)};

    if (p[0] == 0)
        return u_x(p);
    else if (p[0] == 1)
        return -u_x(p);
    else if (p[1] == 0)
        return u_y(p);
    else if (p[1] == 1)
        return -u_y(p);
    else
        return dot(partial, normalize(p - c));
}

const string main_dir = "./main/";

void test_JsonWriter();

bool parseBVP(std::istream &is, BVP2d &solver);

double Conv_Rate(double *norm, int size);

inline void prtNorm(std::ostream &os, double *v);

void SolveBVP(const string &filename);

int main(int argc, char **argv)
{
    for (int i = 1; i < argc; ++i)
        SolveBVP(argv[i]);
    return 0;
}

bool parseBVP(std::istream &is, BVP2d &solver)
{
    Json::Value root;
    Json::Reader reader;
    if (reader.parse(is, root))
    {
        Json::Value center_arr = root["disk"]["center"];
        Point2d center{center_arr[0].asDouble(), center_arr[1].asDouble()};
        double radius = root["disk"]["radius"].asDouble();
        BoundaryCondition BC;
        BC.a = root["BoundaryCondition"]["alpha"].asDouble();
        BC.b = root["BoundaryCondition"]["beta"].asDouble();

        string type = root["BoundaryCondition"]["type"].asString();
        cout << "Solving BVP in " << ((radius > 0) ? "irregular" : "regular") << " domain with boundary condition : " << type << endl;

        solver = BVP2d{
            center, radius, nullptr,
            BC};
        return true;
    }
    else
    {
        // cout << "Invalid input stream!" << endl;
    }
    return false;
}

void SolveBVP(const string &filename)
{
    BVP2d solver;

    cout << SECTOR << endl;
    cout << "Solving the problem specified by : " << filename << endl;

    string prefix = filename.substr(0, filename.size() - 5);
    const string out_dir("./out/" + prefix + "/");
    fs::create_directories("./out/" + prefix);

    ifstream test_file(main_dir + filename);
    if (test_file)
    {
        parseBVP(test_file, solver);
        double alpha = solver.getBC().a;
        double beta = solver.getBC().b;

        BVP2d::Disk D = solver.getD();

        solver.getBC().g = [=](Point2d p)
        { return alpha * u(p) + beta * g_partial(p, D); };

        solver.assignF(f);

        double norm1[4], norm2[4], norm_inf[4];

        for (int i = 0, grid_size = 8; grid_size <= 64; ++i, grid_size *= 2)
        {
            ofstream data(out_dir + to_string(grid_size) + ".txt");
            vector<double> solution = solver.Solve(grid_size);

            for (auto p : solver.getGrid())
            {
                data << p[0] << ' ';
            }
            data << endl;
            for (auto p : solver.getGrid())
            {
                data << p[1] << ' ';
            }
            data << endl;

            vector<double> r;

            int grid_num = solver.getGrid().size();

            double c = solution[0] - u(solver.getGrid()[0]); // for neumann condition

            for (int j = 0; j < solution.size(); ++j)
            {
                r.push_back(solution[j] - u(solver.getGrid()[j]) - c);

                data << solution[j] - c << ' ';

                norm1[i] += fabs(r[j]);
                norm2[i] += r[j] * r[j];
                norm_inf[i] = std::max(norm_inf[i], fabs(r[j]));
            }
            norm1[i] /= grid_num;
            norm2[i] = sqrt(norm2[i] / grid_num);

            data << endl;
            std::cout << "n : " << grid_size << std::endl;
            std::cout << "1-norm : " << norm1[i] << std::endl;
            std::cout << "2-norm : " << norm2[i] << std::endl;
            std::cout << "infinity-norm : " << norm_inf[i] << std::endl;
            std::cout << std::endl;
        }
        ofstream norm(out_dir + "norm.txt");

        prtNorm(norm, norm1);
        prtNorm(norm, norm2);
        prtNorm(norm, norm_inf);

        cout << "Convergence rate: " << endl;
        cout << "1-norm : " << Conv_Rate(norm1, 4) << endl;
        cout << "2-norm : " << Conv_Rate(norm2, 4) << endl;
        cout << "infinity-norm : " << Conv_Rate(norm_inf, 4) << endl;
    }
    else
        cout << "Failed to open the json file!" << endl;
}

double Conv_Rate(double *norm, int size)
{
    double cr = 0;
    for (int i = 1; i < size; ++i)
    {
        cr += (log(norm[i - 1]) - log(norm[i])) / log(2);
    }
    return cr / (size - 1);
}

inline void prtNorm(std::ostream &os, double *v)
{
    for (int i = 0; i < 4; ++i)
        os << v[i] << ' ';
    os << endl;
}

void test_JsonWriter()
{
    Json::Value root, disk, BC;
    disk["center"][0] = (double)0.0;
    disk["center"][1] = (double)0.0;
    disk["radius"] = (double)0.4;
    BC["alpha"] = (double)1.0;
    BC["beta"] = (double)0.0;
    BC["type"] = "Dirichlet";
    root["disk"] = disk;
    root["BoundaryCondition"] = BC;

    Json::StyledWriter writer;
    ofstream out(main_dir + "test.json");
    out << writer.write(root);
    out.close();
}